/*

An example of using the library for EA optimization, using a dummy evaluation function for testing purpose

Author	: Dudy Lim	
Email	: dudy0001@ntu.edu.sg

*/

#include<iostream.h>
#include<string.h>
#include<fstream.h>
#include<vector.h>
#include<math.h>
#include<iomanip.h>
#include "molCluster.h"
#include "population.h"


using std::string;


// please modify this part properly
int clusterSize = 4;	// e.g. (H2O)4 -> 4 molecules (or clusterSize in this context).
int atomPerMolecule = 3;	// e.g. H2O -> 3 atoms.


double extractGauss( ) {
        int i;
        char gaussOutFile[20];
        sprintf(gaussOutFile, "gaussianOut.dat");
        char argu[100];
        sprintf(argu, "/bin/grep \"E\(\" %s > grepOut.dat", gaussOutFile);

        system(argu);

        char dummy[20];
        ifstream grepFile("grepOut.dat");
        for(i=0; i<5; i++)
                grepFile >> dummy;
        grepFile.close();

        double energy;
        return energy = atof( dummy );

}

vector<double> evaluateGaussian( vector<atom>& atomsArg  ) {
        int length =   (clusterSize*atomPerMolecule) ;
        int i, j;

        ofstream gaussianIn("gaussianIn.dat",ios::trunc);
        gaussianIn << "\%chk=W4lm.chk \n\%mem=128MW \n\%nproc=1 \n# b3lyp/6-31+G*\n\n\n0 1 " << endl;

        for(i=0; i<length; i++) {
                gaussianIn << atomsArg[i].label << " ";
                gaussianIn << setprecision(10) << scientific << atomsArg[ i ].x << " ";
                gaussianIn << setprecision(10) << scientific << atomsArg[ i ].y << " ";
                gaussianIn << setprecision(10) << scientific << atomsArg[ i ].z << " ";
                gaussianIn << endl;
        }
        gaussianIn << endl;
        gaussianIn.close();

        system("/opt/gaussian/g03/g03 <  gaussianIn.dat >  gaussianOut.dat");

	vector<double> properties(2);
        properties[0] = extractGauss( );

        return properties;
}



//void evaluateGaussianWithOpt(molCluster& x) {
vector<double> evaluateGaussianWithOpt(  vector<atom>& atomsArg ) {


        int length =   (clusterSize*atomPerMolecule) ;

        int i, j;

        ofstream gaussianIn("gaussianIn.dat",ios::trunc);
        gaussianIn << "\%chk=W4lm.chk \n\%mem=128MW \n\%nproc=1 \n# opt b3lyp/6-31+G*\n\nTitle\n\n0 1" << endl;


        for(i=0; i<length; i++) {
                gaussianIn << atomsArg[ i ].label << " ";
                gaussianIn << setprecision(10) << scientific << atomsArg[ i ].x << " ";
                gaussianIn << setprecision(10) << scientific << atomsArg[ i ].y << " ";
                gaussianIn << setprecision(10) << scientific << atomsArg[ i ].z << " ";
                gaussianIn << endl;
        }
        gaussianIn << endl;
        gaussianIn.close();

        char dummy[256];

        system("/opt/gaussian/g03/g03 <  gaussianIn.dat >  gaussianOpt.dat");

        // get the last structure in the output file, together with the energy and rms, modify molCluster x
        system("~/extractGau-2.0/extractGau -g gaussianOpt.dat -o gaussianOpt.xyz");

        int nStruct;

        //get linenumber of gaussianOpt.xyz
        unsigned long   ulLineCount;
        ifstream gauss("gaussianOpt.xyz");
        ulLineCount = 0;
        while (gauss.getline(dummy, sizeof(dummy), '\n')){
                ulLineCount++;
        }
        gauss.close();
        nStruct = (int)( ceil(ulLineCount/(length+2)) ) - 1;

	
	vector<double> properties(2);

        ifstream gaussianOptXYZ("gaussianOpt.xyz");
        for(i=0; i<nStruct; i++) {
                gaussianOptXYZ >> dummy >> dummy >> properties[0] >> properties[1] >> dummy >> dummy;
                for(j=0; j<length; j++) {
                        gaussianOptXYZ >> atomsArg[ j ].label >> atomsArg[ j ].x >> atomsArg[ j ].y >> atomsArg[ j ].z;
                }
        }

	return properties;
}


//
// A real evaluation function
//
vector<double> myEvaluator( vector<atom>& atomsArg, int optimize ) {


        if( optimize ) return evaluateGaussianWithOpt(  atomsArg );

        else return evaluateGaussian( atomsArg );
}


int main( int arg, char** argv ) {

	srand(time(NULL));
	int i, j;
	int iteration = 20;
	int nElitist = 1;
	int mutationStrength = 2;

	string populationFile("population.dat");
	string molSettingFile("molSetting.dat");

	population parent( populationFile, molSettingFile, myEvaluator );
	population offspring( 1 );
	
	/*	
	cout << "first test" << endl;
	cout << "before evaluate: " << endl;
	for(i=0; i<12;i++) {
                cout <<  parent[0].atoms[i].label << " " << parent[0].atoms[i].x << " " << parent[0].atoms[i].y << "  " << parent[0].atoms[i].z << endl;
        }
	cout << parent[0].properties[0] << " " << parent[0].properties[1] << endl;
	parent[0].evaluate(0);
	cout << "after evaluate: " << endl;
        for(i=0; i<12;i++) {
                cout <<  parent[0].atoms[i].label << " " << parent[0].atoms[i].x << " " << parent[0].atoms[i].y << "  " << parent[0].atoms[i].z << endl;
        }
        cout << parent[0].properties[0] << " " << parent[0].properties[1] << endl;
	

	cout << "\nsecond test" << endl;
	cout << "before evaluate: " << endl;
	for(i=0; i<12;i++) {
                cout <<  parent[0].atoms[i].label << " " << parent[0].atoms[i].x << " " << parent[0].atoms[i].y << "  " << parent[0].atoms[i].z << endl;
        }
	cout << parent[0].properties[0] << " " << parent[0].properties[1] << endl;
        parent[0].evaluate(1);
	cout << "after evaluate: " << endl;
        for(i=0; i<12;i++) {
		cout <<  parent[0].atoms[i].label << " " << parent[0].atoms[i].x << " " << parent[0].atoms[i].y << "  " << parent[0].atoms[i].z << endl;
	}
	cout << parent[0].properties[0] << " " << parent[0].properties[1] << endl;

	*/

	//
	// example of ga+local search
	//

	for(j=0; j<parent.size(); j++) {
		parent[j].evaluate( 0 );
	}	

	for(i=0; i<iteration; i++) {
			
		cout << "Iteration " << i+1 << ", best fitness is " << parent.bestFitness() << endl;
		
		for(j=0; j<offspring.size(); j++) {
			molCluster dad, mum; 
			parent.selectOneRandom( dad );	//note that other types of selections are also available (see population.h)
			parent.selectOneRandom( mum );	
			molCluster& child = offspring[j];
			
			dad.crossover( mum, child );
			child.mutate( mutationStrength );

			child.evaluate( 1 );
		}
		
		parent.replaceMuLambda( 9, offspring );
		
	}

	cout <<	"bestfitness " << parent.bestFitness() << endl;
	cout <<	"worstfitness " << parent.worstFitness() << endl;
	cout <<	"medianfitness " << parent.medianFitness() << endl;
	cout <<	"meanfitness " << parent.meanFitness() << endl;
	cout <<	"stdevfitness "  << parent.stdevFitness() << endl;


	return 0;

}
